// SMBKC Basic three-stage catch-survey-analysis (CSA) model 
// Prerecruit 1: 90-105 mm CL; Prerecurit 2: 105-119 mm CL; Legal: 120+ mm CL

// Constructed by Bill Gaeuman March 2011

// Data for objective function
//  1) trawl survey composition and total abundance (and CV)
//  2) pot survey composition and total abundance (and CV)
//  3) fishery retained average catch weight
//     (fishery total retained weight considered known)
//  4) crab observer composition data from observed count proportions
//
// Quadratic penalties on log-recruit and log-fishing mortality deviations
//
// Directed fishery assumed to occur as pulse at midpoint of season.
//
// Abundances in 1000s of crab (crab per 1000 pot lifts for pot survey estimate).
// Biomasses in 1000s of lb (lb per 1000 pot lifts for pot survey estimate).
// Average weights in lb.
//
// Avoid biomasses except as specifically needed for computing management quantities.
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

DATA_SECTION
 init_int styr            // Beginning year, e.g. 1978
 init_int endyr           // Beginning year, e.g. 2011
 init_int nyrs            // Model time frame in years, e.g. 33
 init_vector wgt(1,3)     // Stage mean weights for any necessary biomass computations
 init_vector hm(1,3)      // Directed and groundfish fixed-gear and trawl fishery handling mortalities

 init_int nyrs_ts                    // Number of years of trawl survey data
 init_ivector yrs_ts(1,nyrs_ts)      // Trawl survey data year indices
 init_matrix ts_data(1,nyrs_ts,1,6)  // Sample size, stage abundance indices, total abundance, CV

 init_int nyrs_ps                     // Number of years of pot survey data
 init_ivector yrs_ps(1,nyrs_ps)       // Pot survey data year indices
 init_matrix ps_data(1,nyrs_ps,1,6)   // Sample size, stage abundance indices, total abundance, CV
 
 init_int nyrs_pf                     // Number of years of directed pot fishery data (other than zero catch)
 init_ivector yrs_pf(1,nyrs_pf)       // Fishery data year indices
 init_matrix pf_data(styr,endyr,1,3)  // Catch wt, avg wt, time to midpoint of fishery 
 
 init_int nyrs_ob                    // Number of years of observer data
 init_ivector yrs_ob(1,nyrs_ob)      // Observer data year indices
 init_matrix ob_data(1,nyrs_ob,1,3)  // Observed stage counts
 //----------------------------------
 
 //Error trap to ensure all data have been read
 init_int eof;
 !! if(eof != 999){cout<<"DATA READING ERROR"<<endl; exit(1);};
 //----------------------------------
 
 ivector yrs(1,nyrs)      // Model years, e.g. 1978, 1979,..., 2010
 
 vector n_ts(1,nyrs_ts)      // Survey and observer data sample sizes [number of male crab >= 90mm CL]
 vector n_ps(1,nyrs_ps)
 vector n_ob(1,nyrs_ob)
 
 vector x_ts(1,nyrs_ts)      // Survey estimated total abundances/biomasses and fishery retained average weight
 vector x_ps(1,nyrs_ps)      // and fishery retained average weight
 vector ret_avg_wgt(styr,endyr)
 
 vector x_ret(styr,endyr)    // Retained catch weight [considered known]
 vector lag_pf(styr,endyr)  // Time to pot fishery [this is zero if no fishery]
 !!lag_pf = 0.0;
 
 vector cv_ts(1,nyrs_ts)    // Survey estimated CVs
 vector cv_ps(1,nyrs_ps)

 matrix p_ts(1,nyrs_ts,1,3)    // Survey and fishery (from observer data) stage proportions
 matrix p_ps(1,nyrs_ps,1,3)    
 matrix p_ob(1,nyrs_ob,1,3) 

 vector sig_ts(1,nyrs_ts);
 vector sig_ps(1,nyrs_ps);
 vector effn_ts(1,nyrs_ts);
 vector effn_ps(1,nyrs_ps);
 vector effn_ob(1,nyrs_ob);
 LOCAL_CALCS
   int k;
   // Vector of years
   yrs.fill_seqadd(styr,1);
   //Extract data
   // Trawl Survey Data
   n_ts  = column(ts_data,1);
   x_ts  = column(ts_data,5);
   cv_ts = column(ts_data,6);
   for(k=1;k<=nyrs_ts;k++) 
     p_ts(k) = --ts_data(k)(2,4)/sum(ts_data(k)(2,4));
  
   // Pot Survey Data
   n_ps  = column(ps_data,1);
   x_ps  = column(ps_data,5);
   cv_ps = column(ps_data,6);
   for(k=1;k<=nyrs_ps;k++) 
     p_ps(k) = --ps_data(k)(2,4)/sum(ps_data(k)(2,4));
  
   // Pot Fishery Data
   for(k=styr;k<=endyr;k++)
   {
     x_ret(k)       = pf_data(k,1);  // retained catch 
     ret_avg_wgt(k) = pf_data(k,2);  // Assumed known  
     lag_pf(k)      = pf_data(k,3);  // This = 0 in years with no fishery
   }
   
   // Observer Data
   n_ob = rowsum(ob_data);
   for(k=1;k<=nyrs_ob;k++) 
     p_ob(k)=ob_data(k)/n_ob(k);
  
   // Trawl and pot-survey abundance data standard deviations for between-years relative weighting
   sig_ts = sqrt( log(square(cv_ts) + 1.0) );
   sig_ps = sqrt( log(square(cv_ps) + 1.0) );
   // Between-years effective sample size relative weighting for composition data--should be moved to data section
   effn_ts = n_ts/max(n_ts);
   effn_ps = n_ps/max(n_ps);
   effn_ob = n_ob/max(n_ob);
 END_CALCS
  !! ad_comm::change_datafile_name("cm.ctl");
  init_int  use_fake        // Flag for whether to use fake data or not...0=not
  init_int  ph_M
  init_int  ph_Q
  init_int  ph_M98
  init_int  ph_Qp
  init_int  ph_logN1o
  init_int  ph_logN2o
  init_int  ph_logN3o
  init_int  ph_logit_p12
  init_int  ph_logit_p23
  init_int  ph_as_ts
  init_int  ph_bs_ts
  init_number  in_as_ts
  init_number  in_bs_ts
  init_int  ph_as_ps
  init_int  ph_bs_ps
  init_int  ph_as_pf
  init_int  ph_bs_pf
  init_int  ph_mean_log_Fpf
  init_int  ph_log_Fpf_dev
  init_int  ph_mean_log_New
  init_int  ph_log_New_dev
  init_vector Lw(1,6)
  init_vector Pw(1,2)
//____________________________________________________________________________
PARAMETER_SECTION
 init_number M(ph_M)    // Natural mortality exept 98/99
 init_number Q(ph_Q)    // Trawl survey catchability
 init_bounded_number M98(0.5,1.5,ph_M98)  // Natural mortality in 98/99
 init_bounded_number Qp(0.1,10.0,ph_Qp)  // Pot survey index proportionality constant
 
 // Log initial stage abundances
 init_bounded_number logN1o(0,15,ph_logN1o)      
 init_bounded_number logN2o(0,15,ph_logN2o)
 init_bounded_number logN3o(0,15,ph_logN3o)
 
 // Logit p12 and p23 transition probabilities
 init_bounded_number logit_p12(-10,10,ph_logit_p12)
 init_bounded_number logit_p23(-10,10,ph_logit_p23)
 
 // Stage 1,2 logistic selectivity function coefficients
 init_bounded_number as_ts(0.01,1.0001,ph_as_ts)      
 init_bounded_number bs_ts(0.01,1.0001,ph_bs_ts)
 init_bounded_number as_ps(0.01,1.0001,ph_as_ps)
 init_bounded_number bs_ps(0.01,1.0001,ph_bs_ps)
 init_bounded_number as_pf(0.01,1.0,ph_as_pf)
 init_bounded_number bs_pf(0.01,1.0,ph_bs_pf)
 
 // Mean log fishing mortality and deviations
 init_number mean_log_Fpf(ph_mean_log_Fpf)    
 init_bounded_dev_vector log_Fpf_dev(styr,endyr,-8.0,8.0,ph_log_Fpf_dev)  
 
 // Mean log recruitment and deviations
 init_bounded_number mean_log_New(5.0,10.0,ph_mean_log_New)        
 init_bounded_dev_vector log_New_dev(styr,endyr,-8.0,8.0,ph_log_New_dev)  
 //-------------------------------------------------------
 
 // Yearly natural mortality [= M98 in year 21 and otherwise = M]
 vector MM(styr,endyr)
 
 // Trawl-survey, pot-survey and directed-fishery stage 1,2 selectivities
 number s1_ts
 number s2_ts
 number s1_ps
 number s2_ps
 number s1_pf // Assumed 0.4 in existing 4-stage model
 number s2_pf // Assumed 0.6 in existing 4-stage model
 
 // Row-stage-to-column-stage transition matrix (molting + growth)
 matrix TM(1,3,1,3)
 
 // Fishing mortalitites [= 0 in years with no fishery]
 vector Fpf(styr,endyr)
 !! Fpf.initialize();
 
 // Model recruitment [note: New(t) contributes to N1(t+1)]
 sdreport_vector New(styr,endyr)
 sdreport_vector Ntot(styr,endyr)
 sdreport_number mmb215
 
 // Yearly stage abundances at beginning of year [survey time]
 matrix N(styr,endyr,1,3)  
 
 // Model predicted fishery stage removal (mortality) numbers [= 0 in years with no fishery]
 matrix R_pf(styr,endyr,1,3)
 !! R_pf.initialize(); // set equal to zero
 
 // Model predicted abundance/biomass measures 
 vector X_ts(1,nyrs_ts)
 vector X_ps(1,nyrs_ps)
 vector Ret_avg_wgt(styr,endyr)
 
 // Model predicted composition measures
 matrix P_ts(1,nyrs_ts,1,3)
 matrix P_ps(1,nyrs_ps,1,3)
 matrix P_ob(1,nyrs_ob,1,3)
  
 // negLoglikelihoods (less additive constants)
 vector L(1,6)
 // negLogPenalties (less additive constants)
 vector P(1,2)
 objective_function_value nll 
//_____________________________________________________________________

INITIALIZATION_SECTION
 // These are also parameter values used to generate fake data
 M  0.18;
 Q  1.0;
 M98 0.5
 Qp 4.0
 logN1o 7.6
 logN2o 7.6
 logN3o 8.0
 logit_p12 1.7
 logit_p23 1.1
 as_ts in_as_ts
 bs_ts in_bs_ts
 as_ps 0.4
 bs_ps 0.6
 as_pf 0.3
 bs_pf 0.8
 mean_log_Fpf -1.5
 mean_log_New 7.6    
//______________________________________________________________________

PRELIMINARY_CALCS_SECTION
 // For validation [comment out to use the .dat file data as is]
  // if (use_fake)
    // mk_fake_data();
//________________________________________________________________________  
   
PROCEDURE_SECTION
 get_numbers(); 
 run_pop_dynamics();
 predict_data();
 calculate_obj_function();
 if (sd_phase())
 {
   for(int i=styr;i<=endyr;i++)
     Ntot(i) = sum(N(i));
   get_mmb215();
 }
 if (mceval_phase())
   write_mcmc();
 // cout<<L<<endl;
   
FUNCTION write_mcmc
  mcout << nll<<" "
        <<M98 <<" "
        <<Qp  <<" "
        <<logN1o <<" "
        <<logN2o <<" "

        <<logN3o <<" "
        <<logit_p12 <<" "
        <<logit_p23 <<" "
        <<as_ts     <<" "
        <<bs_ts     <<" "
        <<as_ps     <<" "
        <<bs_ps     <<" "
        <<as_pf     <<" "
        <<bs_pf     <<" "
        <<New       <<endl;
//________________________________________________________________________
 
FUNCTION get_numbers
 int j;
 
 // Natural mortality for years 1 to nyrs-1
 MM = M; MM(1998) = M98;
 
 // Selectivites from two-parameter logistic curves
 s1_ts = as_ts;
 s2_ts = bs_ts;
 s1_ps = as_ps;
 s2_ps = bs_ps;
 s1_pf = as_pf;
 s2_pf = bs_pf;
 // s1_ts = 1.0/( 1.0+mfexp(-as_ts*(97.0-bs_ts)) );
 // s2_ts = 1.0/( 1.0+mfexp(-as_ts*(112.0-bs_ts)) );
 // s1_ps = 1.0/( 1.0+mfexp(-as_ps*(97.0-bs_ps)) );
 // s2_ps = 1.0/( 1.0+mfexp(-as_ps*(112.0-bs_ps)) );
 // s1_pf = 1.0/( 1.0+mfexp(-as_pf*(97.0-bs_pf)) );
 // s2_pf = 1.0/( 1.0+mfexp(-as_pf*(112.0-bs_pf)) );
 
 //Transition matrix depends on 2 estimated parameters logit_p12, logit_p23 
 dvariable p12, p23;
 p12     = 1.0/( 1.0+mfexp(-logit_p12) );
 p23     = 1.0/( 1.0+mfexp(-logit_p23) );
 TM(1,1) = 1.0-p12; TM(1,2)=p12;     TM(1,3)=0.0;
 TM(2,1) = 0.0;     TM(2,2)=1.0-p23; TM(2,3)=p23;
 TM(3,1) = 0.0;     TM(3,2)=0.0;     TM(3,3)=1.0;
 
 // Estimated fishing mortalities [= 0 in years with no fishery]
  Fpf = mfexp(mean_log_Fpf + log_Fpf_dev); // vector = vector (same shape)

 // Estimated model recruitment [New(t) contributes to N(t+1,1)]
  New = mfexp(mean_log_New + log_New_dev);

 // Initial stage abundances
 N(styr,1)=mfexp(logN1o); N(styr,2)=mfexp(logN2o); N(styr,3)=mfexp(logN3o);
 //__________________________________________________________________________

FUNCTION run_pop_dynamics
 int t;
 dvariable NN1, NN2, NN3;
 dvariable S,D,PS;
 
 for(t=styr;t<=endyr;t++)
 {
   // Survival to directed pot fishery, full-selection fishery mortality fraction, 
   // post-fishery survival
   S  = mfexp(-lag_pf(t)*MM(t)); 
   D  = (1.0-mfexp(-Fpf(t)));
   PS = mfexp(-(1.0-lag_pf(t))*MM(t));
   
   // Calculate fishery removals 
   R_pf(t,1) = N(t,1)*S*D*s1_pf*hm(1);
   R_pf(t,2) = N(t,2)*S*D*s2_pf*hm(1);
   R_pf(t,3) = N(t,3)*S*D;
   
   // Calculate end-of-year pre molt/growth abundances
   NN1 = (N(t,1)*S-R_pf(t,1))*PS;
   NN2 = (N(t,2)*S-R_pf(t,2))*PS;
   NN3 = (N(t,3)*S-R_pf(t,3))*PS;
   
   // Calculate next year's abundances
   if (t<endyr)
   {
     N(t+1,1) = TM(1,1)*NN1+New(t);
     N(t+1,2) = TM(1,2)*NN1+TM(2,2)*NN2;
     N(t+1,3) = TM(2,3)*NN2+NN3;
   }
 }

//________________________________________________________________________

FUNCTION predict_data
 int j;
 // Predicted average retained weight is (assumed known) retained biomass/number retained
 for(j=styr;j<=endyr;j++)
   Ret_avg_wgt(j) = x_ret(j)/R_pf(j,3); // total catch biomass over predicted N...

 // Predicted trawl survey total abundance and proportions
 for(j=1;j<=nyrs_ts;j++)
 {
   X_ts(j)   = N(yrs_ts(j),1)*s1_ts + N(yrs_ts(j),2)*s2_ts + N(yrs_ts(j),3);
   P_ts(j,1) = N(yrs_ts(j),1)*s1_ts/X_ts(j);
   P_ts(j,2) = N(yrs_ts(j),2)*s2_ts/X_ts(j);
   P_ts(j,3) = N(yrs_ts(j),3)/X_ts(j);
 }
 X_ts = Q*X_ts;
 
 // Predicted pot-survey total abundance and proportions
 for(j=1;j<=nyrs_ps;j++)
 {
   X_ps(j) = N(yrs_ps(j),1)*s1_ps + N(yrs_ps(j),2)*s2_ps + N(yrs_ps(j),3);
   P_ps(j,1)= N(yrs_ps(j),1)*s1_ps/X_ps(j);
   P_ps(j,2)= N(yrs_ps(j),2)*s2_ps/X_ps(j);
   P_ps(j,3)= N(yrs_ps(j),3)/X_ps(j);
 }
 X_ps = X_ps/Qp;

 // Predicted observer proportions using stage removals [after accounting for handling mortality]
 for(j=1;j<=nyrs_ob;j++)
 {
   P_ob(j,1) = R_pf(yrs_ob(j),1)/hm(1);
   P_ob(j,2) = R_pf(yrs_ob(j),2)/hm(1);
   P_ob(j,3) = R_pf(yrs_ob(j),3);
   P_ob(j) = P_ob(j) / sum(P_ob(j));
 }
//___________________________________________________________________________

FUNCTION calculate_obj_function
 int j;
           
 // 1. Retained avg weight normally distributed about predicted retained avg weight
 L(1) = 0.5*norm2(ret_avg_wgt-Ret_avg_wgt);

 // 2. Trawl suvey abundance lognormally distributed about predicted value
 L(2) = 0.5*norm2(elem_div(log(x_ts)-log(X_ts),sig_ts));
 
 // 3. Pot survey abundance lognormally distributed about predicted value
 L(3) = 0.5*norm2(elem_div(log(x_ps)-log(X_ps),sig_ps));

 // 4. Trawl survey proportions are multinomial wrt predicted proportions
 L(4) = -effn_ts*rowsum(elem_prod(p_ts,log(P_ts + 0.01)));  
 
 // 5. Pot survey proportions are multinomial wrt predicted proportions
 L(5) = -effn_ps*rowsum(elem_prod(p_ps,log(P_ps + 0.01)));
 
 // 6. Observer proportions are multinomial wrt predicted proportions
 L(6) = -effn_ob*rowsum(elem_prod(p_ob,log(P_ob + 0.01)));

 // 1. Model recruit deviations implies CV on R as about 0.4...
 P(1) = 0.5*norm2(log_New_dev);

 // 2. Directed fishery log fishing mortality deviations
 P(2) = 0.5*norm2(log_Fpf_dev);

 // Objective function
 L = elem_prod(L,Lw); // apply weights here so finals are known (in report)
 P = elem_prod(P,Pw); // apply weights here so finals are known (in report)
 nll =  sum(P) + sum(L);
 if (!last_phase())
   nll += norm2(Fpf-.1)*12.5; // a modest penalty to keep it real during earlier estimation phases.

REPORT_SECTION
 report<<"F"<<endl;
 report<<Fpf<<endl;
 report<<"Initial numbers"<<endl;
 report<<N(styr)<<endl;
 report<<"TS, PS, and fishery selectivities"<<endl;
 report<<s1_ts<<"  "<<s2_ts<<endl;
 report<<s1_ps<<"  "<<s2_ps<<endl;
 report<<s1_pf<<"  "<<s2_pf<<endl;
 report<<"Transition Matrix"<<endl;
 report<<TM<<endl;
 report<<"PS Q"<<endl;
 report<<Qp<<endl;
 report<<"1998 mortality"<<endl;
 report<<M98<<endl;
 report<<"Predicted and reported average weights"<<endl;
 report<<Ret_avg_wgt<<endl;
 report<<ret_avg_wgt<<endl;
 report<<"Recruits"<<endl;
 report<<New<<endl;
 write_R();

FUNCTION write_R
 R_report<<"$Fpf"<<endl;
 R_report<<Fpf<<endl;
 R_report<<"$init_N"<<endl;
 R_report<<N(styr)<<endl;
 dvector tmp(1,20);
 tmp.fill_seqadd(0,8); // fill with 0,8, ... 160 mm
// set up the initial concentrations of the two reactants fo
 R_report<<"$mm"<<endl;
 R_report<<tmp<<endl;
 R_report<<"$sel_ts"<<endl;
 R_report<<as_ts<<" "<<bs_ts<<" 1 "<<endl;
 R_report<<"$sel_ps"<<endl;
 R_report<<as_ps<<" "<<bs_ps<<" 1 "<<endl;
 R_report<<"$sel_pf"<<endl;
 R_report<<as_pf<<" "<<bs_pf<<" 1 "<<endl;
 R_report<<"$trans_mat"<<endl;
 R_report<<TM<<endl;
 R_report<<"$q_ps"<<endl;
 R_report<<Qp<<endl;
 R_report<<"$M98"<<endl;
 R_report<<M98<<endl;
 R_report<<"$pred_aw"<<endl;
 R_report<<Ret_avg_wgt<<endl;
 R_report<<"$obs_aw"<<endl;
 R_report<<ret_avg_wgt<<endl;
 R_report<<"$recs"<<endl;
 R_report<<New<<endl;
 R_report<<"$yrs_pf"<<endl;
 R_report<<yrs_pf<<endl;
 R_report<<"$yrs_ts"<<endl;
 R_report<<yrs_ts<<endl;
 R_report<<"$x_ts"<<endl;
 R_report<<x_ts<<endl;
 R_report<<"$X_ts"<<endl;
 R_report<<X_ts<<endl;
 R_report<<"$yrs_ps"<<endl;
 R_report<<yrs_ps<<endl;
 R_report<<"$x_ps"<<endl;
 R_report<<x_ps<<endl;
 R_report<<"$X_ps"<<endl;
 R_report<<X_ps<<endl;
 R_report<<"$p_ps"<<endl;
 R_report<<p_ts<<endl;
 R_report<<"$P_ps"<<endl;
 R_report<<P_ts<<endl;
 R_report<<"$p_ts"<<endl;
 R_report<<p_ts<<endl;
 R_report<<"$P_ts"<<endl;
 R_report<<P_ts<<endl;
 R_report<<"$n_ts"<<endl;
 R_report<<n_ts<<endl;
 R_report<<"$n_ts"<<endl;
 R_report<<n_ts<<endl;
 R_report<<"$n_ps"<<endl;
 R_report<<n_ps<<endl;
 R_report<<"$n_ob"<<endl;
 R_report<<n_ob<<endl;
 R_report<<"$lag_pf"<<endl;
 R_report<<lag_pf<<endl;
 R_report<<"$L"<<endl;
 R_report<<L<<endl;
 R_report<<"$P"<<endl;
 R_report<<P<<endl;
 R_report<<"$N"<<endl;
 R_report<<N<<endl;
//____________________________________________________________________________

FUNCTION get_ADFG_inputs
 // This function outputs legal and mature male abundance and MMB at survey time 
 // for use in ADF&G harvest strategy.
 
 cout<<"Legal Abundance "<<N(nyrs,3)<<endl;
 cout<<"Mature Male Abundance "<<sum(N(nyrs)(2,3))<<endl;
 cout<<"Mature Male Biomass "<<N(nyrs)(2,3)*wgt(2,3)<<endl;
 //____________________________________________________________________________

FUNCTION get_MSST
// This function outputs a baseline MSY based on estimates of MMB215 for 
// determination of MSST = 1/2MSY.
//______________________________________________________________________________

FUNCTION get_mmb215
// This function outputs estimate of mmb215 in most recent fishery year for
// determination against MSST of overfished status. Must declare as sdreport_number
// in main code if want estimated standard error or as likeprof_number if want that.
 
 dvariable n2, n3;
 
 n2 = ((N(endyr,2)*mfexp(-lag_pf(endyr)*MM(endyr))-R_pf(endyr,2))*mfexp(-(0.63-lag_pf(endyr))*MM(endyr)));
 n3 = ((N(endyr,3)*mfexp(-lag_pf(endyr)*MM(endyr))-R_pf(endyr,3))*mfexp(-(0.63-lag_pf(endyr))*MM(endyr)));
 mmb215 = n2*wgt(2)+n3*wgt(3);
 // cout<<mmb215<<endl;
//______________________________________________________________________________

FUNCTION get_OFL
// This function outputs OFL and F_OFL using F35% and fixed-point methodology.


//_______________________________________________________________________________

GLOBALS_SECTION
 #include <math.h>
 #include <admodel.h>
 ofstream R_report("cm_R.rep");
 ofstream mcout("cm_mcmc.rep");
//___________________________________________________________________________



//________________________________________________________________________

FUNCTION dvector multi3(const int& seed, const int& size, const dvector& p)
// Returns 3-vector of multinomial counts given seed, sample size, and
// 3-vector p, which is first normalized to produce probabilities. 
 {
   p /= sum(p);
   random_number_generator rng(seed);
   dvector u(1,size);
   double a, b;
   dvector n(1,3);
   int j;
   n = 0.0; a = p(1); b = p(1) + p(2);
   u.fill_randu(rng);
   for(j=1;j<=size;j++)
   {
     if(u(j)<a) n(1)+=1.0; else
     if(u(j)<b) n(2)+=1.0; else
     n(3)+=1.0;
   }
   return(n);
 }

FUNCTION mk_fake_data
// Overwrites input data with simulated data using following initial parameter values:
 // M = 0.18
 // Q = 1.0
 // M98 = 1.0
 // Qp = 4.0
 // logN1o = logN2o = 7.601; logN3o = 8.001
 // logit_p12 = 1.735; logit_p23 = 1.099
 // as_ts = 0.02 ; bs_ts = 40 ; as_ps = 0.05; bs_ps = 60; as_pf = 0.1; bs_pf = 100;
 // mean_log_Fpf = -1.5
 // mean_log_New = 7.601 
 // log Fpf and log New deviations drawn from normal with standard deviations 0.3 and 0.4
 
 int j,t, seed;
 seed = 67593;
 random_number_generator rng(seed);
 double X;
 int jj,nn;
 dvector pp(1,3);
 
 // Derived parameters needed for simulation
 MM=0.18; MM(21)=1.0;
 
 N(styr,1)=2000.0; N(styr,2)=2000.0; N(styr,3)=3000.0;
 
 double p12=0.85;
 double p23=0.75;
 TM(1,1)=1.0-p12; TM(1,2)=p12;     TM(1,3)=0.0;
 TM(2,1)=0.0;     TM(2,2)=1.0-p23; TM(2,3)=p23;
 TM(3,1)=0.0;     TM(3,2)=0.0;     TM(3,3)=1.0;
 
 s1_ts = 1.0/( 1.0+mfexp(-0.02*(97.0-40.0)) );
 s2_ts = 1.0/( 1.0+mfexp(-0.02*(112.0-40.0)) );
 s1_ps = 1.0/( 1.0+mfexp(-0.05*(97.0-60.0)) );
 s2_ps = 1.0/( 1.0+mfexp(-0.05*(112.0-60.0)) );
 s1_pf = 1.0/( 1.0+mfexp(-0.1*(97.0-100.0)) );
 s2_pf = 1.0/( 1.0+mfexp(-0.1*(112.0-100.0)) );
 
 // Directed fishery fishing mortalities [= 0 for non-fishery years]
 dvector Fdevs(1,nyrs_pf);
 Fdevs.fill_randn(rng);
 for(j=1;j<=nyrs_pf;j++)
   Fpf(yrs_pf(j))=mfexp(-1.5 + Fdevs(j)*0.3);   

 // Model recruits
 dvector Rdevs(styr,endyr);
 Rdevs.fill_randn(rng); 
 // for(j=styr;j<nyrs;j++)
   // New(j) = mfexp(7.601+Rdevs(j)*0.4);
 New = mfexp(7.601+Rdevs*0.4);
 
 double S,D,PS,NN1,NN2,NN3;
 for(t=styr;t<endyr;t++) // only to yr before last...
 {
    // Survival to directed pot fishery, full-selection fishery mortality fraction, 
    // post-fishery survival
    S  = value(mfexp(-lag_pf(t)*MM(t)));
    D  = value((1.0-mfexp(-Fpf(t))));
    PS = value(mfexp(-(1.0-lag_pf(t))*MM(t)));
    
    // Calculate fishery removals 
    R_pf(t,1) = N(t,1)*S*D*s1_pf*hm(1);
    R_pf(t,2) = N(t,2)*S*D*s2_pf*hm(1);
    R_pf(t,3) = N(t,3)*S*D;
    
    // Calculate end-of-year pre molt/growth abundances
    NN1 = value((N(t,1)*S-R_pf(t,1))*PS);
    NN2 = value((N(t,2)*S-R_pf(t,2))*PS);
    NN3 = value((N(t,3)*S-R_pf(t,3))*PS);
    
    // Calculate next year's abundances
    N(t+1,1) = TM(1,1)*NN1+New(t);
    N(t+1,2) = TM(1,2)*NN1+TM(2,2)*NN2;
    N(t+1,3) = TM(2,3)*NN2+NN3;
 }
 
 // Known retained catch biomass and new ret_avg_wgt data from normal draw around old data value
 dvector wgt_devs(styr,endyr);
 wgt_devs.fill_randn(rng);
 for(j=styr;j<=endyr;j++)
 {
   x_ret(j) = value( R_pf(j,3)*ret_avg_wgt(j) );
   ret_avg_wgt(j) += wgt_devs(j)*0.05;
 }

 // Trawl survey composition and total abundance data
 dvector epsilon_ts(1,nyrs_ts);
 epsilon_ts.fill_randn(rng);
 for(j=1;j<=nyrs_ts;j++)
 {
   jj = yrs_ts(j);
   nn = (int)n_ts(j);
   pp(1) = value(s1_ts*N(jj,1));
   pp(2) = value(s2_ts*N(jj,2));
   pp(3) = value(N(jj,3));
   pp = pp/sum(pp);
   p_ts(j) = multi3(seed+1,nn,pp);
   p_ts(j) = p_ts(j)/sum(p_ts(j));

   X = value(Q*(s1_ts*N(jj,1)+s2_ts*N(jj,2)+N(jj,3)));
   epsilon_ts(j) *= sqrt(log(square(cv_ts(j))+1.0));
   x_ts(j) = X*mfexp(epsilon_ts(j));  
 }

 // Pot survey composition and total abundance data
 dvector epsilon_ps(1,nyrs_ps);
 epsilon_ps.fill_randn(rng);
 for(j=1;j<=nyrs_ps;j++)
 {
   jj = yrs_ps(j);
   nn = (int)n_ps(j);
   pp(1) = value(s1_ps*N(jj,1));
   pp(2) = value(s2_ps*N(jj,2));
   pp(3) = value(N(jj,3));
   pp = pp/sum(pp);
   p_ps(j) = multi3(seed+2,nn,pp);
   p_ps(j) = p_ps(j)/sum(p_ps(j));
   
   X = value((s1_ps*N(jj,1)+s2_ps*N(jj,2)+N(jj,3))/Qp);
   epsilon_ps(j) *= sqrt(log(square(cv_ps(j))+1.0));
   x_ps(j) = X*mfexp(epsilon_ps(j));       
 }
 
 // Observer data proportions using stage removals [sample size n_ob is sum of stage counts]
 for(j=1;j<=nyrs_ob;j++)
 {
   nn = (int)n_ob(j);
   pp(1) = value( R_pf(yrs_ob(j),1)/hm(1) );
   pp(2) = value( R_pf(yrs_ob(j),2)/hm(1) );
   pp(3) = value( R_pf(yrs_ob(j),3) );
   pp = pp/sum(pp);
   p_ob(j) = multi3(seed,nn,pp);
   p_ob(j) = p_ob(j)/sum(p_ob(j));
 }
//____________________________________________________________________________


